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ABSTRACT 

The Parker or field line tangling model of coronal heating is studied comprehensively via long- 
time high-resolution simulations of the dynamics of a coronal loop in cartesian geometry within the 
framework of reduced magnetohydrodynamics (RMHD). Slow photospheric motions induce a Poynting 
flux which saturates by driving an anisotropic turbulent cascade dominated by magnetic energy. In 
physical space this corresponds to a magnetic topology where magnetic field lines are barely entangled, 
nevertheless current sheets (corresponding to the original tangential discontinuities hypothesized by 
Parker) are continuously formed and dissipated. 

Current sheets are the result of the nonlinear cascade that transfers energy from the scale of con- 
vective motions (~ 1,000 km) down to the dissipative scales, where it is finally converted to heat 
and/or particle acceleration. Current sheets constitute the dissipative structure of the system, and 
the associated magnetic reconnection gives rise to impulsive "bursty" heating events at the small 
scales. This picture is consistent with the slender loops observed by state-of-the-art (E)UV and X-ray 
imagers which, although apparently quiescent, shine bright in these wavelengths with little evidence 
of entangled features. 

The different regimes of weak and strong MHD turbulence that develop, and their influence on 
coronal heating scalings, are shown to depend on the loop parameters, and this dependence is quanti- 
tatively characterized: weak turbulence regimes and steeper spectra occur in stronger loop fields and 
lead to larger heating rates than in weak field regions. 

Subject headings: MHD — Sun: corona — Sun: magnetic fields — turbulence 



1. INTRODUCTION 

In a previous letter l|Rappazzo et al.l |2007( ) we de- 
scribed simulations, within the framework of RMHD in 
cartesian geometry, aimed at solving the Parker field-lin e 
tangling (coronal heating) problem (Par k"erl[l972t fl988). 
We also developed a phenomenological model for nonlin- 
ear interactions, taking into account the inertial photo- 
spheric line-tying effect, which explained how the aver- 
age coronal heating rate would depend on the only free 
parameter present in the simulations, namely the ratio 
of the coronal loop Alfven crossing time and the pho- 
tospheric eddy turnover time. This paper is devoted to 
a more detailed discussion of the numerical simulations 
and of the relationship between this work, the original 
Parker conjecture, and the nanoflare scenario of coronal 
heating. 

Parker's book (|Parker|[l994f ) is devoted to an exami- 
nation of the basic theorem of magnetostatics, namely 
that the lowest available energy state of a magnetic field 
in an infinitely conducting fluid contains surfaces of tan- 
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gential discontinuity, or current sheets. It is Parker's 
conjecture that the continuous footpoint displacement of 
coronal magnetic field lines must lead to the development 
of such discontinuities as the field continuously tries to 
relax to its equilibrium state, and it is the dynamical 
interplay of energy accumulation via footpoint motion 
and the bursty dissipation in the forming current sheets 
which gives rise to the phenomenon of the high temper- 
ature solar corona, heated by the individual bursts of 
reconnection, or nanoflares. 

What then does turbulence have to do with the 
nanoflare heating scenario? Parker himself strongly criti- 
cizes the use of the "t" word, the formation of the current 
sheets being due in his opinion to the "requirement for 
ultimate static balance of the Maxwell stresses". But 
what better way is there to describe the nonlinear global 
dynamics of a magnetically dominated plasma in which 
the formation of an equilibrium state containing current 
sheets is the inevitable asymptotic state (once the pho- 
tospheric driver is turned off)? 

The striving of the global magnetic field toward a state 
containing current sheets must occur through local vio- 
lations of the force-free condition, the induction of local 
flows, the collapse of the currents into ever thinner layers: 
a nonlinear process generating ever smaller scales. From 
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the spectral point of view, a power law distribution of 
energy as a function of scale is expected, even though 
the kinetic energy is much smaller than the magnetic en- 
ergy. The last two statements are clear indications that 
the word turbulence provides a correct description of the 
dynamical process. 

A final important issue is whether the overall dissi- 
pated power tends to a finite value as the resistivity and 
viscosity of the coronal plasma become arbitrarily small. 
That this must be the case is easy to understand (see 
§ 13. 3[) . For suppose that for an arbitrary, continuous, 
foot-point displacement the coronal field were only to 
map the foot-point motion, and that there were no non- 
linear interactions, i.e. the Lorentz force and convective 
derivatives were negligible everywhere. In this case, the 
magnetic field and the currents in the corona would then 
grow linearly in time, until the coronal dissipation at 
the scale of photospheric motions balanced the forcing. 
The amplitudes of the coronal fields and currents would 
then be inversely proportional to resistivity (eqs. (|30[) - 
(|31|) ). and the dissipated power, product of resistivity 
and square of the current, would also scale as the in- 
verse power of the resistivity (eq. ([33]) ). In other words, 
the smaller the resistivity in the corona, the higher the 
power dissipated would be. But the amplitudes can not 
become arbitrarily large, because non-linear effects inter- 
vene to stop the increase in field amplitudes, increasing 
the effective dissipation at a given resistivity. Since the 
power can not continue to increase monotonically as the 
resistivity is decreased, it is clear that at some point non- 
linear interactions must limit the dissipated power to a 
finite value, regardless of the value of the resistivity. Fi- 
nite dissipation at arbitrarily small values of dissipative 
coefficients is another definition of a turbulent system. 

All this assuming that a statistically stationary state 
may be reached in a finite time, a question closely re- 
lated to the presence of finite time singularities in 3D 
magnetohydrodynamics. It now appears that magnetic 
field relaxation in an unforced situation does not lead to 
the development of infinitely thin current sheets in a fi- 
nite time, but rather the cu rrent development appears t o 
be only exponential in time (|Grauer and Marliani |[2000h ■ 
In forced numerical simulations, as the ones we will de- 
scribe in detail here, this is a moot point: for all intents 
and purposes a statistically stationary state is achieved 
at a finite time independent of resistivity for sufficiently 
high resolution. In fact, even if the growth is expo- 
nential, we can estimate that the width of the current 
sheets reaches the meter-scale in a few tens Alfven cross- 
ing times T4. A typical value is T4 = 40 s, so that this 
initial time is not only finite, but also short compared 
with an active region timescale. Once the steady state 
has been reached this phenomenon is no longer impor- 
tant. The nonlinear regime is in fact characterized by the 
presence of numerous current sheets, so that while some 
of them are being dissipated others are being formed, 
and a statistical steady state is maintained. 

It therefore seems that the Parker field-line tan- 
gling scenario of coronal heating may be described 
as a particular instance of magnetically dominated 
MHD turbulence. Numerous analytical and numerical 
models of this process have been presented in the 
past, each discussing in some detail a spects of the 
general problem as presented above (|Parkerl I1972L 
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The numerical simulations presented here bring closure 
to the original question as posed in cartesian geometry 
by Parker, starting from a uniform axial magnetic field 
straddling from one boundary plane to another, subject 
to continuous independent footpoint motions in either 
photosphere. This does not imply that we have fully 
solved the coronal heating problem as due to footpoint 
dragging by the photospheric velocity field. 

A number of relevant effects have been neglected: first, 
the field line expansion between the photosphere and 
corona, which, if the photospheric flux is confined to 
bundles in granular and supergranular network lanes, 
would allow the mapping of the photospheric velocity 
field to the coronal volume to contain discontinuities. 
We are presently carrying out a dedicated set of sim- 
ulations to capture this effect. Second, the projection of 
the 3D photospheric velocity to 2D coronal base motions 
parallel to the photosphere also introduces compressibil- 
ity in the forcing flow, again neglected here. Third, we 
have considered stationary photospheric flows. The ef- 
fect of a finite eddy-turnover time in the flow w a s cons id- 
ered in lEinaudi et alJ <|1996f ): IGeorgoulis et all (|1998f l in 
2 dimensions, and in the "3 d i mensio nal" shell model 
calculations of iBuchlin fc Vellil (|2007ft . These showed 
that time-dependence does not change things substan- 
tially provided the flow pattern does not contain degen- 
erate symmetries, a fact confirmed by shorter simula- 
tions we defer to a future paper. Finally, we do not ad- 
dress the more realistic case of a single photosphere with 
curved cor onal loops, such as the simul ations presented 
recently by I Gudiksen fc Nordlundl (|2005D . While this ap- 
proach has advantages when investigating the coronal 
loop dynamics within its coronal neighborhood, modeling 
a larger part of the solar corona numerically drastically 
reduces the number of points occupied by the coronal 
loops. At the moment the very low resolution attainable 
with this kind of simulations does not allow the develop- 
ment of turbulence with a well-developed inertial range. 
The transfer of energy from the scale of convection cells 
~ 1000 km toward smaller scales is inhibited, because the 
smaller scales are not resolved (their linear resolution is 
~ 500 km). Thus, these simulations have not been able 
to shed light on the detailed coronal statistical response 
nor on the different regimes which may develop and how 
they depend on the coronal magnetic field crossing time 
and the photospheric eddy turnover time. 

In § [2] we introduce the coronal loop model, whose 
properties are qualitatively analyzed in § [3] The results 
of our simulations are described in § [4j and their turbu- 
lence properties are analyze in more detail in § [5l Finally 
in §[6] we summarize and discuss our results. 

2. PHYSICAL MODEL 

A coronal loop is a closed magnetic structure threaded 
by a strong axial field, with the footpoints rooted in 
the photosphere. This makes it a strongly anisotropic 
system, as measured by the relative magnitude of the 
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Alfven velocity associated with the axial magnetic field 
va ~ 2000 kms -1 compared to the typical photosphcric 
velocity u p h ~ 1 kms -1 . 

We study the loop dynamics in a simplified Carte- 
sian geometry, neglecting field line curvature, i.e. the 
toroidality of loops. Our loop is a "straightened out" 
box, with an orthogonal square cross section of size I 
(along which the x-y directions lie), and an axial length 
L (along the z direction) embedded in an axial homoge- 
neous uniform magnetic field Bo = fio e 2 . This simpli- 
fied geometry allows us to perform simulations with both 
high numerical resolution and long-time duration. 

In § 12.11 we introduce the equations used to model the 
dynamics, while in § 12.21 we give the boundary and initial 
conditions used in our numerical simulations. 

2.1. Governing Equations 

The dynamics of a plasma embedded in a strong ax- 
ial magnetic field are well described by the equations 
of reduced MHD (RMHD) jlKadomtsev '& Poguts"i ll974t 
IStrausslll976tlMontgomervlll982fl . 

These equations are valid for a plasma with small ratio 
of kinetic to magnetic pressures, in the limit of a large 
loop-aspect ratio (e = l/L <C 1, L being the length of the 
loop and I being the minor radius of the loop) and of a 
small ratio of poloidal to axial magnetic field (b±/Bo < 
e). In dimensionless form they can be written as: 



du± 
~dt 



Ren 



(1) 



db±_ 



(b ± -V ± )u ± -{u ± -V_ 



du i 



(-1) 



n+1 



dz 



Re r , 



V2n i 



0. 



0. 



(3) 



where u_l and h± are the components of the velocity and 
magnetic fields perpendicular to the mean field, and p is 
the kinetic pressure. The gradient operator likewise has 
only components in the x-y plane perpendicular to the 
axial direction z, i.e. 



Vi =e : 



d_ 

dx 



d_ 

dy' 



(4) 



while the dynamics in the planes is coupled to the axial 
direction through the linear terms cx d z . 

To render the equation non dimensional magnetic fields 
have first been expressed in velocity units by dividing by 
^/Airpo (where po is a density supposed homogeneous and 
constant), i.e. considering the associated Alfven veloci- 
ties (b — ► b/^/Anpo), and then both velocity and magnetic 
fields have been normalized to a typical photospheric ve- 
locity Uph] lengths and times have been expressed in units 
of the perpendicular length of the computational box t 
and its related "eddy turnover time" t± = £/u p h- 

As a result, in equations (H])-©, the linear terms cx 
d z are multiplied by the dimensionless Alfven velocity 
C A = V A/ U ph, i-e the ratio between the Alfven velocity 



associated with the axial magnetic field va = -Bo/ V^Po? 
and the photospheric velocity u p h- 

The index n is called dissipativity: the diffusive terms 
adopted in equations dU)- ([2j) correspond to ordinary dif- 
fusion for n = 1 and to so-called hyperdiffusion for n > 1 . 
When n — 1 the V\/ Re diffusive operator is recovered, 
so that Re\ — Re — Re m corresponds to the kinetic 
and magnetic Reynolds number (considered of equal and 
uniform value): 



Re = 



Po ^ u 



ph 



Re m = 



ph 



(5) 



where viscosity v and resistivity r\ are taken to be con- 
stant and uniform (c is the speed of light) . 

We have performed numerical simulations with both 
n = 1 and n — 4. Hyperdiffusion is used because with a 
limited resolution the diffusive timescales associated with 
ordinary diffusion are small enough to affect the large 
scale dynamics and render very difficult the resolution of 
an inertial range, even with a grid with 512x512 points in 
the x-y plane (the highest resolution grid we used for the 
plane). The diffusive time r„ at the scale A associated 
with the dissipative terms used in ([I])-© is given by: 



Re n A 2 



(G) 



While for n = 1 the diffusive time decreases relatively 
slowly towards smaller scales, for n — 4 it decreases 
far more rapidly. This allows to have longer diffu- 
sive timescales at large spatial scales and similar dif- 
fusive timescales at the resolution scale. Numerically 
we require that the diffusion time at the resolution scale 
A m in = where N is the number of grid points, to 

be of the same order of magnitude for both normal and 
hyperdiffusion, i.e. 



Rei 



Re n 
TV 2 ™ 



Re n ~ Re\ AT 2 (" _1 ) (7) 



For instance a numerical grid with N = 512 points which 
requires a Reynolds number Rei — 800 with ordinary 
diffusion, can implement Re^ ~ 10 19 , removing diffusive 
effects at the large scales, and allowing (if present) the 
resolution of an inertial range. 

The numerical integration of the RMHD equations (fT])- 
([3]) is substantially simplified by using the potentials of 
the velocity (ip) and magnetic field (ip), 



u ± =Vx(tpe z ), b ± = V x (ipe z ) , 



(8) 



linked to vorticity and current by u = — V^c/j and j = 

-vi</>. 

We solve numericall y equations ([Tl)-(l3l) writ ten in terms 
of the potentials (see iRappazzo etalT ( 2007f) ) in Fourier 
space, i.e. we advance the Fourier components in the x-y 
directions of the scalar potentials ip and -0- Along the 
z direction no Fourier transform is performed so that 
we can impose non-periodic boundary conditions (speci- 
fied in § 12. 2[) . and a central second-order finite difference 
scheme is used. In the x-y plane a Fourier pseudospec- 
tral method is implemented. Time is discretized with a 
third-order Runge-Kutta method. 

We use a computational box with an aspect ratio of 
10, which spans 



0<x,y<l, 



< z < 10. 
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2.2. Boundary and Initial Conditions 

As boundary conditions at the photosphcric surfaces 
(z = 0, L) we impose two independent velocity pat- 
terns, intended to mimic photospheric motions, made up 
of large spatial scale projected convection cell flow pat- 
terns constant in time. The velocity potential at each 
boundary is given by: 



1 



E 



sin 



27rVfc 2 + I 2 
{kx + ly) + 2n £_ ki 



(10) 



We excite all the wave number values (fc, I) € 1? in- 

1/2 

eluded in the range 3 < (fc 2 + I 2 ) < 4, so that the 
resulting average injection wavenumber is k c ~ 3.4, and 
the average injection scale £ c , the convection cell scale, 
is given by £ c — £/k c . a k i and are two sets of ran- 
dom numbers whose values range between and 1, and 
are independently chosen for the two boundary surfaces. 
The normalization adopted in eq. (|10[) sets the value of 

the corresponding velocity rms (see eq. (JSJ) to 1/V2, i.e. 



\ \ dxdy (i 
Jo Jo 



(11) 



At time t = no perturbation is imposed inside the 
computational box, i.e. h± = u± = 0, and only the axial 
magnetic field Bq is present: the subsequent dynamics 
are then the effect of the photospheric forcing (flTJ|) on 
the system, as described in the following sections. 

3. ANALYSIS 

In order to clarify aspects of the linear and nonlinear 
properties of the RMHD system, we provide an equiv- 
alent form of the equations (HJ-©. In terms of the 
Elsasser variables z ± = u± ± bj_, which bring out the 
basic symmetry of the equations in terms of parallel and 
anti-parallel propagating Alfven waves, they can be writ- 
ten as 



dz+ 
dt 



~~dt 



■ V, 



CA 



dz+ 



i-iy 



Re„ 



v±p, 



(12) 



dz- 



(-1) 



n+1 



Re,, 



ViF, (13) 
(14) 

where P = p + h\/2 is the total pressure, and is linked 
to the nonlinear terms by incompressibility (|14l) : 



VlP = -J2{d, 

i,3=l 



(15) 



In terms of the Elsasser variables z = uj_ ± b 

0,L 



velocity pattern at upper or lower boundary surface 

becomes the constraint z + + z = 2u°^ L at that bound- 
ary. Since, in terms of characteristics (which in this case 



are simply z,^ themselves), we can specify only the in- 
coming wave (while the outgoing wave is determined by 
the dynamics inside the computational box) , this velocity 
pattern implies a reflecting condition at the top (z = L) 
and bottom (z = 0) planes: 



2u 



2u 



at z = 0, 



at z = L. 



(16) 
(17) 



The linear terms (<x d z ) in equations (fT2)) - (fT3|) give rise 
to two distinct wave equations for the fields, which 
describe Alfven waves propagating along the axial direc- 
tion z. This wave propagation, which is present dur- 
ing both the linear and nonlinear regimes, is responsible 
for the continuous energy influx on large perpendicular 
scales (see eq. (TTU)) ) from the boundaries into the loop. 
The nonlinear terms (z T • V^)z ± are then responsible 
for the transport of this energy from the large scales to- 
ward the small scales, where energy is finally dissipated, 
i.e. converted to heat and/or particle acceleration. 

A well-known important feature of the nonlinear terms 
in equations p2 |l - (fT4|l is the absence of self-coupling, i.e. 
only counterpropagating waves interact non-linearly, and 
if one of the two fields z ± is zero, there are no non- 
linear interactions at all. This fact, i.e. that counter- 
propagating wave-packets may interact only while they 
are crossing e ach other, lies at the basis of the s o-called 
Alfven effect ( Iroshnikovl 1 1964t iKraichnanl 1 1965T ) . which 
ultimately renders the nonlinear timcscales longer and 
slows down the dynamics. 

From this description three different timescales arise 
naturally: ta, T p h and r„;. r_4 = L/va is the crossing 
time of the Alfven waves along the axial direction z, i.e. 
the time it takes for an Alfven wave to cover the loop 
length L. r p h ~ 5 m is the characteristic time associ- 
ated with photospheric motions, while r„; is the nonlin- 
ear timescale. 

For a typical coronal loop T4 <^ r p h ■ For instance for a 
coronal loop long L = 40, 000 km and with an Alfven ve- 
locity vj, = 2, 000 fcms -1 we obtain ta — 20 s, which is 
small compared to r p h ~ 5 m = 300 s. This is the reason 
we carried out simulations with a photospheric forcing 
constant in time (see eq. (|10jl). i.e. for which formally 

T ph = OO. 

In the RMHD ordering the nonlinear timescale r„; is 
bigger than the Alfven crossing time ta- As we shall see 
this ordering is maintained during our simulations and 
we will give analytical estimates of the value of r„; as a 
function of the characteristic parameters of the system. 

An important feature of equations (fT2|) -(|14 p that we 
will use to generalize our results is that, apart from the 
Reynolds numbers, there is only one fundamental non- 
dimensional parameter: 



/ = 



4 VA 

Luph 



(18) 



Hence all the physical quantities which result from the 
dynamical evolution, e.g. energy, Poynting flux, heating 
rate, timescales, etc., must depend on this single param- 
eter /. 

3.1. Energy Equation 
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From equations ([l])-©, with n = 1, and considering 
the Reynolds numbers equal, the following energy equa- 
tion can be derived: 

|Qui + ibi)=-V.S-±(jW), ,19) 

where S = B x (u x B) is the Poynting vector. As ex- 
pected the energy balance of the system described by 
cq. (fT9|) is due to the competition between the energy 
(Poynting) flux flowing into the computational box and 
the ohmic and viscous dissipation. Integrating eq (I19|) 
over the whole box the only relevant component of the 
Poynting vector is the component along the axial direc- 
tion z, because in the x-y plane periodic boundary con- 
ditions are used and their contribution to the Poynting 
flux is null. As B — c A e z + b± and u = , this is given 

by 

S z = S ■ e z = -c A (u ± ■ b±) . (20) 

Considering that the velocity fields at the photospheric 
boundaries are given by Uj_ and Uj_, for the integrated 
energy flux we obtain 

S = c A I da (ui • b ± ) -c A I da « • bj_) . (21) 

Jz=L Jz=0 

The injected energy flux therefore depends not only on 
the photospheric forcing and the axial Alfven velocity 
(which have fixed values), but also on the value of the 
magnetic fields at the boundaries, which is determined 
by the dynamics of the system inside the computational 
box: the injection of energy depends on the nonlinear 
dynamics which develops, and viceversa. 

The simplified topology investigated in this paper, 
i.e. a strong axial magnetic field whose footpoints are 
dragged by 2D orthogonal motions applies to regions 
where emerging flux may be neglected. Consider the ax- 
ial component of the velocity u z field carrying new mag- 
netic field (bjH into the corona. The associated Poynt- 
ing flux is 

Sf=(h e /)\ z . (22) 

This flux component is negligible when S e J < S Z: i.e., 
since all the components of the photospheric velocity 
fields are of the same order, u z ~ u p h, when 

(b e l) 2 <B bT b - (23) 

In § 15.21 we give an estimate of the value of the field 
jjturb g enera ted by the field-line dragging, and will be 

able to quantify for which value of b e / the emerging flux 
can be neglected. 

3.2. Linear Stage 

For t < r n i nonlinear terms can be neglected. Ne- 
glecting also the diffusion terms, which play no role on 
large scales, equations ([Tj)-(|3]) reduce to two simple wave 
equations. Coupled with the boundary conditions (|10[) 
the solution for times longer than the crossing time r A 
reads: 

b±(x,y,z,t) = [u L (x,2/)-u (x,2/)] — , (24) 

T A 

ui(x, y, z, t) = u L (x, y)j + u°(x, y) (l - J) . (25) 



This shows that A) the loop velocity field is bounded 
by the imposed photospheric fields and b) the magnetic 
field grows linearly in time, uniform along the loop, while 
mapping the photospheric velocity field in the perpen- 
dicular planes. Therefore, for a generic set of velocities 
u L and u°, the resulting magnetic fields (f2"4"|) -([25" f give 
rise to non-vanishing forces in the perpendicular planes 
which grow quadratically in time, becoming dynamicall y 
important after a certain interval ( Buchlin & Vclli 2007) . 

There exists a (singular) set of velocity forcing pat- 
terns, for which the generated coronal field has a van- 
ishing Lorentz force. For simplicity consider u L = 0: 
in terms of potentials it follows that ip = —(p°t/T A and 
tp = (fi° (1 — z/L) (where = V x (ip° e 2 )). In this 
case both and uj^ are proportional to Vj_ x (ip° e z ). 
The condition for the vanishing of nonlinear terms then 
becomes 

V (V 2 ip°) x Vip° = 0, with tp° = tp° (x, y) . (26) 

This condition is then satisfied by those fields for which 
the laplacian is constant along the streamlines of the 
field. As u> = — V 2 Lp this is equivalent to the statement 
that the vorticity is constant along the streamlines. This 
condition is in general not verified, unless very symmet- 
ric functions are chosen, e.g. in cartesian geometry by 
any ID function like ip° — f(x), and in polar coordinates 
by any radial function (p° — g(r). 

Generally speaking even in such peculiar configurations 
non-linear interactions will arise due to the onset of insta- 
bilities. We defer discussion of these extreme examples 
to a subsequent paper, the random photospheric fields 
(I10p discussed here always giving rise to non-vanishing 
forces . 

Inserting the linear evolution fields (|24p in the expres- 
sion for the integrated energy flux (f2Tj) . we find 

S = c A [da\u L -u°\ 2 ■ — , (27) 

J T A 

i.e. the Poynting flux S grows linearly in time until such 
a time that non-linear interactions set in. 

A similar l inear analysis was already performed by 
Parker (1988), who noted that if this is the mechanism 
responsible for coronal heating, then the energy flux 
S z ~ S/£ 2 must approach the value S z ~ 10 7 erg cm 2 s" 1 
necessary to sustain an active region before a saturat- 
ing mechanism, magnetic reconnection of singular cur- 
rent sheets in Parker's language, takes over. 

In fact however the value reached by S z depends on 
the nonlinear dynamics, its value self- consistently deter- 
mined by solving the nonlinear problem. An S z too small 
compared with observational constraints would then rule 
out the Parker model. 

3.3. Effects of Diffusion 

The linear solution ([24|) - (|25|l has been obtained with- 
out taking into account the diffusive terms. This is jus- 
tified, given the large value of the Reynolds numbers for 
the solar corona. But numerically it can be important. 
At very low resolution diffusion is so important that little 
or no nonlinear dynamics develop and the system reaches 
a balance between the photospheric forcing and diffusion 
of the large scale fields. 

One can attempt to bypass the non-linear problem 
by adopting a much smaller "turbulent" value of the 
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Fig. 1. 

boundary forcing at the bottom plane z 
In lighter vortices the velocity field is directed anti- 
clockwise while in darker vortices it is directed clock- 
wise. The cross-section shown in the figure is roughly 
4000 x 4000 km 2 , where the typical scale of a convective 
cell is 1000 km. 



Reynolds number ijHevvaerts fcPriest I fl99 2). For this 
"ad hoc" value of the Reynolds number the average dissi- 
pation would be the same as in the high Reynolds number 
active turbulence limit. Linearizing equation ([2]) (with 
n = 1 and Re\ = Re), we obtain 



<9bj 
dt 



= CA- 



dz 



(28) 



Taking into account that the forcing velocities are dom- 
inated by components at the injection scale i c (see 

eq. (HO])), the relation V^p = — (2tt/£ c ) ip, where £ c = 
£/k c with the average wavenumber k c ~ 3.4, is approx- 
imately valid. Integrating then eq. (|2"5|) over z and di- 
viding by the length L, we obtain for averaged along 
z: 



db±_ 
~~dt 



L 



[u L (x,y)-u° (x,y)] 



£ 2 Re 



(29) 



Indicating with w. p h = u L — u°, with Tn — £ 2 Re/ (2ir) 2 
the diffusive time-scale and with T4 = L/ca the Alfven 
crossing time, the solution is given by: 



b± {x,y,t) = u ph (x,y) — 



1 — exp 



, (30) 



\j (•'■■//• = U-ph (•''•//) ( 



1 — exp 



(31) 



So that the magnetic energy Em and the ohmic dissipa- 




Fig. 2. — : Streamlines of the velocity field u^, the 
boundary forcing at the top plane z = L for run A. The 
numerical grid has 512x512 points in the x-y planes, with 
a linear resolution of ~ 8 km. 



tion rate J are given by 



1 



E M = - I d 6 x b 



1 & t 2 ( T n 
-t Lu ph I — 



TA 



1 — exp 



J = — 



1 

Re 



d 3 a; f = 



£ 2 Lu 2 ph 



1 — exp 



t 



t 

t-r 



(32) 



(33) 



where u p h is the rms of u p h, and with the rms of the 
boundary velocities u° and u L fixed to 1/2 pT|) we have 
u p h ~ 1 • Both total magnetic energy and ohmic dis- 
sipation (|33p grow quadratically in time for time smaller 
than the resistive time Tn, while on the diffusive time 
scale they saturate to the values 



1/6 rO. 



jTisat 
b M = 



u* vh R<? 



L (27rfc c ) 



J 3 



u P h Re 



L (27rfc c ) 



(34) 



written explicitly in terms of the loop parameters and 
Reynolds number. 

Magnetic energy saturates to a value proportional to 
the square of both the Reynolds number and the Alfven 
velocity, while the heating rate saturates to a value that 
is proportional to the Reynolds number and the square 
of the axial Alfven velocity. We have also used equa- 
tions (f3"2"l) - (l3"3"l) as a check in our numerical simulations, 
and during the linear stage, before nonlinearity sets in 
they are well satisfied. 

From equation (|3"2"]) - (f3"3")) we can estimate the satura- 
tion time as the time at which the functions (l32l)-J33l) 
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Run 




Tlx X fly X Tlx 




Re, Rei 


+ /'T A 
f'VUCLX / ' J\ 


A 


200 


512 x 512 x 200 


1 


8 • 10 2 


548 


13 


200 


256 x 256 x 100 


1 


4 • 10 2 


1061 


C 


200 


128 x 128 x 100 


1 


2 ■ 10 2 


2172 


D 


200 


128 x 128 x 100 


1 


1 • 10 2 


658 


E 


200 


128 x 128 x 100 


1 


1 • 10 1 


1272 


F 


50 


512 x 512 x 200 


4 


3 ■ 10 20 


196 


G 


200 


512 x 512 x 200 


4 


10 19 


453 


H 


400 


512 x 512 x 200 


4 


10 20 


77 


I 


1000 


512 x 512 x 200 


4 


10 19 


502 



TABLE 1: Summary of the simulations, ca is the ax- 
ial Alfven velocity and n x x n y x n z is number of points 
for the numerical grid, n is the dissipativity, n = 1 indi- 
cates normal diffusion, n — 4 hyperdiffusion. Re (= Re\) 
or Re4 indicates respectively the value of the Reynolds 
number or of the hyperdiffusion coefficient (see eq. (|f 2p ~ 
(jf 3p ). The duration of the simulation t max jTj\ is given 
in Alfven crossing time unit ta — L/va- 

reach 2/3 of the saturation values. It is approximately 
given by 

2t 2 Re 



2tk = 



(2irk c 



(35) 



In the next section we describe the results of our sim- 
ulations, which investigate the linear and nonlinear dy- 
namics. The average values may be used in conjunc- 
tion with[33]to define the equivalent turbulence Reynolds 
number. 

4. NUMERICAL SIMULATIONS 

In this section we present a series of numerical simula- 
tions, summarized in Table [1] modeling a coronal layer 
driven by a forcing velocity pattern constant in time. On 
the bottom and top planes we impose two independent 
velocity forcings as described in § 12.21 which result from 
the linear combination of large-scale eddies with random 
amplitudes, normalized so that the rms of the photo- 
spheric velocity is u p h ~ 1 kms . For each simulation a 
different set of random amplitudes is chosen, correspond- 
ing to different patterns of the forcing velocities. A re- 
alization of this forcing with a specific choice (run A) of 
the random amplitudes is shown in Figures [T][2] 

The length of a coronal section is taken as the unitary 
length. As we excite all the wavenumbers between 3 and 
4, and the typical convection cell scale is ~ 1, 000 km, this 
implies that each side of our section is roughly 4, 000 km 
long. Our typical grid for the cross-sections has 512x512 
grid points, corresponding to ~ I28 2 points per convec- 
tive cell, and hence a linear resolution of ~ 8 km. 

Between the top and bottom plates a uniform magnetic 
field B — Bq e z is present. The subsequent evolution is 
due to the shuffling of the footpoints of the magnetic field 
lines by the photospheric forcing. 

In the different numerical simulations, keeping fixed 
the cross-section length (~ 4, 000 km) and axial length 
(~ 40, 000 km), we explore the behavior of the system for 
different values of ca, i.e. the ratio between the Alfven 
velocity associated with the axial magnetic field and the 
rms of the photospheric motions (density is supposed uni- 
form and constant). 

Nevertheless, as shown in (fT5)l the fundamental param- 
eter is / = £ c va/ Liiph, so that changing ca — VA/u p h 
is equivalent to explore the behaviour of the system for 



different values of /, where the same value of / can be re- 
alized with a different choice of the quantities, provided 
that the RMHD approximation is valid, i.e. we are de- 
scribing a slender loop threaded by a strong magnetic 
field. 

We also perform simulations with different numerical 
resolutions, i.e. different Reynolds numbers, and both 
normal [n = 1) and hyper-diffusion (n = 4). 

The qualitative behaviour of the system is the same 
for all the simulations performed. In the next section 
wc describe these qualitative features in detail for run A, 
and then describe the quantitative differences found in 
the other simulations. 

4.1. Run A 

In this section we present the results of a simula- 
tion performed with a numerical grid with 512x512x200 
points, normal (n = 1) diffusion with a Reynolds number 
Re = 800, and the Alfven velocity va = 200 km s^ 1 cor- 
responding to a ratio ca — VA/uph = 200. The stream- 
lines of the forcing velocities applied in the top (z = L) 
and bottom (z = 0) planes are shown in Figures[T][2l The 
total duration is roughly 550 axial Alfven crossing times 
(t a = L/v A ). 

Plots of the total magnetic and kinetic energies 

E M = \JdVb 2 ± , E K = ^JdVu 2 ± , (36) 
and of the total ohmic and viscous dissipation rates 



J 



n = ±-ldvtu 2 , c$7) 



along with the incoming energy rate (integrated Poynting 
flux) S (see eq. (|2ip). are shown in Figures [3]|4] At the 
beginning the system has a linear behavior (see eqs. ([24]) - 
(|2"B|) , and (f2"T)) ). characterized by a linear growth in time 
for the magnetic energy, the Poynting flux and the elec- 
tric current, which implies a quadratic growth for the 
ohmic dissipation cx (i/r^) 2 , until time t ~ 6r_4, when 
nonlinearity sets in. We can identify this time as the 
nonlinear timescale, i.e. r„; ~ 6T4. The timescales of 
the system will be analyzed in more details in £15.51 

After this time, in the fully nonlinear stage, a statisti- 
cally steady state is reached, in which the Poynting flux, 
i.e. the energy that is entering the system for unitary 
time, balances on time average the total dissipation rate 
(J + f2). As a result there is no average accumulation of 
energy in the box, beyond what has been accumulated 
during the linear stage, while a detailed examination of 
the dissipation time series (see inset in Figure shows 
that the Poynting flux and total dissipations are decor- 
related around dissipation peaks. 

In the diffusive case from eqs. (j3"2")l - (1331) . with the values 
of this simulation we would obtain r sat ~ 50 ta, ~ 
6100 and J sat ~ 7100; all values well beyond those of the 
simulation. A value of Re = 85 would fit the simulated 
average dissipation, while Re = 140 would approximately 
fit the average magnetic energy. In any case this would 
only fit the curves, but the physical phenomena would 
be completely different, as we describe in the following 
sections. 

An important characteristic of the system is the mag- 
netic predominance for both energy and dissipation (Fig- 
ures [3] and 2]) . In the linear stage (§ I3.2[) while the mag- 
netic field grows linearly in time, the velocity field does 
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Fig. 3. — : Run A: High-resolution simulation with 
v A /u ph = 200, 512x512x200 grid points and Re = 800. 
Magnetic (Em) and kinetic (Ek) energies as a function 
of time (t_4 = L/v A is the axial Alfven crossing time). 

not, and its value is roughly the sum of the boundary 
forcing fields. The physical interpretation is that be- 
cause we are bending the axial magnetic field with a con- 
stant forcing, as a result the perpendicular magnetic field 
grows linearly in time, while the velocity remains limited. 
More formally this is a consequence of the fact that, while 
on the perpendicular magnetic field no boundary condi- 
tion is imposed, the velocity field must approach the im- 
posed boundary values at the photosphere both during 
the linear and nonlinear stages. 

In Figure [5] the 2D averages in the x-y planes of the 
magnetic and velocity fields and of the ohmic dissipation 
j 2 /Re, are plotted as a function of z at different times. 
These macroscopic quantities are smooth and present al- 
most no structure along the axial direction. The reason is 
that every disturbance or gradient along the axial direc- 
tion, at least considering the large perpendicular scales 
(for the small scales behavior see § [U ) , is smoothed out 
by the fast propagation of Alfven waves along the axial 
direction, their propagation time T4 is in fact the fastest 
timescale present (in particular T4 < t„j), and then the 
system tends to be homogeneous along this direction. 

The predominance of the ohmic over the viscous dissi- 
pation is due to the fact that, as we show in the next sec- 
tions, the dissipative structures are current sheets, where 
magnetic reconnection takes place. 

The phenomenology described in this section is general 
and we have found it in all the simulations that we have 
performed, in particular we have always found that in 
the nonlinear stage a statistical steady state is reached 
where energies fluctuate around a mean value and to- 
tal dissipation and Poynting flux on the average balance 
while on small timescales decorrelate. In particular, to 
check the temporal stability of these features, which are 
fully confirmed, we have performed a numerical simula- 
tion (run C) with the same parameters as run A, but 
with a lower resolution (128x128x100), a Reynolds num- 
ber Re — 200 and a longer duration (t ~ 2, 000 t a ). On 
the opposite the average levels of the energies and of total 
dissipation depend on the parameters used as we describe 
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Fig. 4. — : Run A: The integrated Poynting flux S dy- 
namically balances the Ohmic (J) and viscous (Q) dissi- 
pation. Inset shows a magnification of total dissipation 
and S for 150 < t/r A < 250. 

in the next sections. 

Before describing these features, in the next section 
we describe the current sheets formation, their temporal 
evolution and other properties. 

4.1.1. Current Sheets, Magnetic Reconnection, Global 
Magnetic Field Topology and Self- Organization 

The nonlinear stage is characterized by the presence of 
current sheets elongated along the axial direction (Fig- 
ures (18a)-(18b)), which exhibit temporal dynamics and 
are the dissipative structures of the system. We now 
show that they are the result of a nonlinear cascade. 
Figure [6] shows the time evolution of the first 11 modes 
of magnetic energy for the first 20 crossing times r A for 
run A. During the linear stage the magnetic field is given 
by eq. and is the mapping of the difference between 
the top (z = 10) and bottom (z = 0) photosphcric ve- 
locities u L (x, y) — u°(x, y), whose streamlines are shown 
in Figure 17a. The field lines of the orthogonal magnetic 
field in the midplane (z = 5) at time t = 0.63 T4 are 
shown in Figure 17b, and as expected they map the ve- 
locity field. The same figure shows in colour the axial 
current j. As shown by eq. (|24|) (taking the curl) the 
large scale motions that we have imposed at the photo- 
sphere induce large scale currents in all the volume and, 
as described in the previous section, if there was not 
a nonlinear dynamics a balance between diffusion and 
forcing would be reached, where no small scale would 
be formed and the magnetic field would always map the 
photosphcric velocities. 

As time proceeds the magnetic field grows and a cas- 
cade transfers energy from the large scales, where the 
photosphcric forcing (|10[) injects energy at the wavenum- 
bers n = 3 and 4, to the small scales (Figured]). In 
physical space this cascade corresponds to the collapse 
of the large scale currents which lead to the formation 
of current sheets, as shown in Figures 17c and 17d. In 
Figures 17e and 17f we show the magnetic field lines at 
time t = 18.47 r_4, in the fully nonlinear stage, with 
respectively the axial component of the current j and 
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Fig. 5. — : Run A: 2D averages in the x-y planes of the 
ohmic dissipation j 2 /Re, the magnetic and velocity fields 
bj_, and u^, as a function of z. The different colours 
represent 10 different times separated by At = 50 in 
the interval 30 < t < 480 t a . 
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Fig. 6. — : Run A: First 11 magnetic energy modes as 
a function of time covering the first 20 Alfven crossing 
times T4. Photosphcric motions inject energy at n = 3 
and 4. 



of the vorticity uj. The resulting magnetic topology is 
quiet complex, X and Y-points are not in fact easily dis- 
tinguished. They are distorted and very often a com- 
ponent of the magnetic field orthogonal to the current 
sheet length is present, so that the sites of reconnection 
are more easily identified by the corresponding vorticity 
quadrupoles. As shown in Figures 17e and 17f, the more 
or less distorted current sheets are always embedded in 
quadrupolar structures for the vorticity, a characteris- 
tic maintained throughout the whole simulation, and a 
clear indication that nonlinear magnetic reconnection is 
taking place. 

Figures 18a and 18b show a view from the side and 
the top of the 3D current sheets at time t = 18.47 ta- 
When looked from the side the current sheets, which are 
elongated along the axial direction, look space filling, but 
the view from the top shows that the filling factor is 
actually small (see also Figure 17). 

Another aspect of the dynamics is self- organization: 
while until time t = 4.79 ta the magnetic field lines are 
still approximately a mapping of the photospheric veloc- 
ities, in the fully nonlinear stage they depart from it and 
have an independent topology that evolves dynamically 
in time (see the associated movie for the time evolution 
covering 40 crossing times from ~ 508 ta up to ~ 548 ta', 
notations and simulation are the same used in Figure 17). 
The reason for which the photospheric forcing does not 
determine the spatial shape of the magnetic field lines 
is due to the bigger value of the rms of the magnetic 
field b± =< b^ > : / 2 in the volume respect to the rms of 
the photosphcric forcings u p h =< (u*[ — u^) 2 > 1 / 2 ^ 1 
(eqs. (16)-C!ID). 

This means that the contribution to the dynamics of 
the Alfvenic perturbations propagating from the bound- 
ary are much smaller, over short periods of time, than the 
self-consistent non-linear evolution due to the magnetic 
fields inside the domain, and therefore can not deter- 
mine the topology of the magnetic field. For run A and 
G, both with ca — 200, the ratio is b^/u p h ~ 6 and it 



increases up to b±/u p h ~ 27 in run I with ca = 1000. On 
the other hand these waves continuously transport from 
the boundaries the energy that sustain the system in a 
magnetically dominated statistically steady state. 

All the facts presented in this section, and the proper- 
ties of the cascade and of the resulting current sheets in 
presence of a magnetic guide field outlined in § El lead to 
the conclusion that the current sheets do not generally 
result directly from a geometrical misalignment of neigh- 
boring magnetic field lines stirred by their footpoints mo- 
tions, but that they are the result of a nonlinear cascade 
in a self- organized system. 

Although the magnetic energy dominates over the ki- 
netic energy, the ratio of the rms of the orthogonal mag- 
netic field over the axial dominant field Bq is quite small. 
For c A = 200, 400 and 1000 it is ~ 3%, so that the aver- 
age inclination of the magnetic fieldlines respect to the 
axial direction is just ~ 2°, it is only for the lower value 
c A = 50 that b±/B ~ 4% and the angle is ~ 4°. The 
field lines of the total magnetic field at time t = 18.47rx 
are shown in Figures 18c and 18d. The computational 
box has been rescaled for an improved viewing, and to at- 
tain the original aspect ratio the box should be stretched 
10 times along the axial direction. The magnetic topol- 
ogy for the total field is quiete simple, as the line ap- 
pear slightly bended. It is only in correspondence of the 
small scale current-sheets that field lines on the opposite 
side may show a relative inclination. But as the current 
sheets are very tiny (and their width decreases at higher 
Reynolds numbers), they occupy only a very small frac- 
tion of the volume, so that the bulk of the magnetic field 
lines appears only slightly bended. 

It is often suggested, or implicitly assumed, that cur- 
rent sheets are formed because the magnetic field line 
footpoints are subject to a random walk. The complexity 
of the footpoint trajectory would then be a necessary in- 
gredient. In fact it would give rise to a complex topology 
for the coronal magnetic field, leading either to tangled 
field lines which would then release energy via fast mag- 
netic reconnection, or to turbulence. So that the "com- 
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plexity" of the footpoint motions would be responsible 
for the "complex" dynamics in the corona. 

On the opposite our simulations show that this system 
in inherently turbulent, and that "simple" footpoint mo- 
tions give rise to turbulent dynamics characterize by the 
presence of an inertial range (§E]) and dynamical current 
sheets. In fact our photospheric forcing velocities (Fig- 
ures [T][2]) are constant in time and have only large-scale 
components (eq. (|10[) ). so that the footpoint motions are 
"ordered" and do not follow any random walk. During 
the linear stage this gives rise to a magnetic field that 
grows linearly in time (eq. (|24| ) and that is a mapping 
of the velocity fields (see eq. ([54)) and Figures 17a and 
17b), i.e. both the magnetic field and the current have 
only large-scale components. The footpoint motions of 
our photospheric velocities never bring two magnetic field 
lines close to one another, i.e. they never geometrically 
produce a current sheet. Current sheets are produced on 
an ideal timescale, the nonlinear timescale, by the cas- 
cade. Furthermore, as we show in the next section, the 
statistically steady state that characterizes the nonlinear 
stage results from the balance at the large-scales between 
the injection of energy and the flow of this energy from 
the large scales toward the small scales, where it is finally 
dissipated. 

As the system is self-organized and the magnetic en- 
ergy increases at higher values of the axial magnetic 
fields, very likely different static or time-dependent (with 
the characteristic photospheric time ~ 300 s) forcing 
functions, will not be able to determine the spatial shape 
of the orthogonal magnetic field. In our more realistic 
simulation with ca — 1000 the ratio b±/u p h is in fact 
~ 27. Other forcing functions are currently being inves- 
tigated, and time-dependent forcing functions are likely 
to modulate with their associated timescale the rms of 
the system, like total energy and dissipation. 

5. TURBULENCE 

Before analyzing in detail further aspects of our simu- 
lations, namely inertial spectra, anisotropics and scaling 
laws, let us briefly justify the statement that the time- 
dependent Parker problem, i.e. the dynamics of a mag- 
netofluid threaded by a strong axial field whose foot- 
points are stirred by a velocity field, is an MHD turbu- 
lence problem. 

The fact that at the large orthogonal scales the Alfvcn 
crossing time ta is the fastest timescale so that dur- 
ing the linear stage the fields evolves as (j2~4")) - ([2"5")) , 
means that the photosphere's role is to contribute an 
anisotropic magnetic forcing function that stirs the fluid, 
with an orthogonal length typical of the convective cells 
(~ 1000 km) and an axial length is given by the loop 
length L. 

Typica l ly, for ced MHD turbulence simulations (e.g. see 
Biskamp ( 2003) and references therein) are performed us- 
ing a 3-periodic numerical cube with a volumetric forcing 
function which mimics some physical process injecting 
energy at the large scales. 

Solutions (f2~4")) - (|25[) can be approximately obtained in- 
troducing themagnetic forcing function F TO in equa- 
tion {2J 

u L (x. iA — u°(x. u) 



and implementing 3-periodic boundary conditions in our 
elongated (0 < x,y < 1, < z < L) computational 
box. During the linear stage this forcing would give rise, 
apart from the small velocity field (|25|) , to the same mag- 
netic field. During the nonlinear stage, as ta < t„i, 
it would still give rise to a similar injection of energy. 
This propert y was the basis for the body of previous 2D 
calculations (lEinaudi et all 119961 : iDmitruk et"afl 119981 : 
iGeorgoulis et al.lll998h ~ 

In particular the photospheric motions imposed at the 
boundaries for the Parker problem take the place of, and 
represent a different physical realizations of the forcing 
function generally used for the 3-periodic MHD turbu- 
lence box. In the Parker model, the equivalent forc- 
ing stirs the magnetic field, whiile in standard simula- 
tions the forcing stirs both velocity and magnetic fields 
or mostly the velocity field. The main differences be- 
tween "standard" MHD turbulence simulations and the 
problem at hand are that a) the peculiarity of the low- 
frequency photospheric forcing leads to magnetic energy 
largely dominating over the kinetic energy in the system 
b) the forcing involves line-tying of the magnetic field 
with 3-periodic boundary conditions. Line-tying inhibit 
the inverse cascade for the magnetic field, as described 
later in this section (§ 15 .4|) . Equivalently, one may say 
that line-tying hinders magnetic reconnection by render- 
ing it less energetically favorable due to the increased 
field line-curvature it requires compared to the unbound 
system. This property is fundamental to the anomalous 
scaling laws and enhanced overall heating rates that will 
be found below. 

In MHD, the cascade takes place preferentially in 
planes orthogo n al to the local mean magnetic field 
(jShebalin et alJ (|1983h ). The small scales formed are 
not uniformly distributed in this plane, rather they are 
organized in dynamical current-vortex sheets extended 
along the direction of the local main field. These cur- 
rent sheets with associated quadrupolar vorticity fila- 
men ts form the dissipative structures of MHD t urbulence 
(e.g. iB'iskamp fc Mulled (|2000f l jBiskanTpl (|2003h and ref- 
erences therein). In our case, because the axial field is 
strong, the current sheets are elongated along the axial 
direction to the point of being quasi-uniform along the 
loop axis (Figure 18). 

5.1. Spectral Properties 

In order to investigate inertial range spectra, we have 
carried out four simulations (runs F, G, H and I in Ta- 
ble [T|) with a resolution of 512x512x200 grid points using 
a mild power (n — 4) for hyperdiffusion (p~2]) - (p~3|) . 

In turbulence the fundamental physical fields are the 
Elsasser variables z ± = ± . Their associated ener- 
gies 

E ± = \JdV (z±) 2 , (39) 

are linked to kinetic and magnetic energies Ek , Em and 
to the cross helicity H c 

H c = ^JdVu ± -b±_ (40) 

by 
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Fig. 7. — : Run G: Ratio between cross-helicity H c and 
total energy E as a function of time. H c <C E shows 
that the system is in a regime of balanced turbulence. 



Nonlinear terms in equations (|12[) - (|15[) are symmetric un- 
der the exchange z + <-> z~ , so as substantially are also 
boundary conditions (fT6|) -(fTT l) . given that the two forc- 
ing velocities are different but have the same rms values 
(= l/\/2). It is then expected that H c « £ so that 
none of the two energies prevails E + ~ E~ ~ E, where 
E = Ek + Em is total energy. In Figure [7] the ratio 
H c JE is shown as a function of time for run G. Cross 
helicity has a maximum value of 5% of total energy, and 
its time average is ~ 1%, and similar values are found for 
all the simulations. Furthermore perpendicular spectra 
of E and E ± in simulations F, G, H and I, overlap each 
other, so that as expected we can also assume that 



5zl 



5zs ~ 6z\, 



(42) 



where 8z\ is the rms value of the Elsasser fields z ± at 
the perpendicular scale A. 

In the following we always consider the spectra in the 
orthogonal plane x-y integrated along the axial direc- 
tion z, unless otherwise noted. Furthermore as they are 
isotropic in the Fourier k x -k y plane, we will consider the 
integrated ID spectra, so that for total energy 



E = - [ dz 
2 Jo 



dz i 



dx dy (it 



E 



u\ 2 



|6| 2 ) (fc,*)=5>n, 

n 

71 = 1,2,... (43) 



where, similarly to eq. (|10|) , n indicates "rings" in k- 
space. Figure [8] shows the total energy spectra E n av- 
eraged in time, obtained from the hyperdiffusive simula- 
tions F, G, H and I with dissipativity n = 4 (eqs. ([T2j) - 
fB]l) and respectively c A = 50, 200, 400 and 1000. An 
inertial range displaying a power law behaviour is clearly 
resolved. The spectra visibly steepens increasing the 
value of c_4, with spectral index ranging from 1.8 for 
= 50 up to ~ 2.7 for ca = 1000. The spectra are 
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Fig. 8. — : Total energy spectra as a function of the 
wavenumber n for simulations F, G, H and I. To higher 
values of ca — v^/uph, the ratio between the Alfven and 
photospheric velocities, correspond steeper spectra, with 
spectral index respectively 1.8, 2, 2.3 and 2.7. 



clearly always steeper than the well known (strong) MHD 

inertial range turbulence spectra k ± 5 ^ or k^ 2 . 

This steepening is certainly not a numerical arti- 
fact: the use of hyperdiffusion gives rise to a hump at 
high wave-number y alues, known as the bottleneck ef- 
fect (jFalkovichl 1 1 994h . which when present flattens the 
spectra. Furthermore we use the same value o f diss i- 
pativity (n = 4) used by iMaron fc Goldreichl (pOOlh . 
who find the same IK spectral slope (—3/2), also con- 
fir med in recent highe r -resol ution simulations performed 
by iMiiller fc Grappinl (|2005t ) with standard n = 1 diffu- 
sion. 

In our simulations, a hump or flattening at high 
wavenumbers is best visible in run H with ca = 400, 
which might be due to the bottleneck effect, but a more 
probable interpretation involves a transition from weak 
to strong turbulence at the smaller scales within the in- 
ertial range, which requires a preliminary discussion of 
strong vs. weal turbulence in MHD. 

Recently a lot of progress has been made in the 
understanding MHD turb ulence both in the condi- 
tion of so-called strong jGoldreich fc Sridhari H995L 
19971: IChofc Vishniad [20001: iBiskamp fc Mulled |200~ 
Muller et all 120031 : IMiiller fc Grappinl 12005: B oldvre 



20051 120061: iMason et alj 120061) and weak turb ulence 



dNg fc Bhattacharied 119971: iGoldreich fc Sridhari H997l 
iGaltier et all 1200a iGaltier fc Chandranl 120061 ). Weak 
turbulence has been investigated mainly through ana- 
lytical methods. The total energy spectrum can be char- 
acterized by a fcj 2 power law, which is easily found phe- 
nomenologically by considering that the Alfven effect oc- 
curs along the field while the cascade proceed s in the 
orthogonal direction (|Ng fc B hattacharice 1997). While 
our MHD simulations, even with our line-tying boundary 
conditions and anomalous energetic regime (b dominat- 
ing over u except at the smallest scales), confirm the 
presence of the A; J spectrum for a range of loop param- 
eters, steeper spectra are also found nearly reaching kj 3 , 
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Fig. 9. — : Run I: Snapshot of the 2D spectrum 
E(nj_,n z ) in bilogarithmic scale at time t ~ 145 r^. 
ni and n z are respectively the orthogonal and axial 
wavenumbers. The 2D spectrum is shown as a function 
of n± and n z + 1, to allow the display of the n 2 = 
component. 



clearly linked to the strength of the axial field Bq an ef- 
fect we discuss more in detail in the following subsection. 

The formation of an inertial range is crucially related 
to the anisotropy of the cascade, where a relationship be- 
tween spectral extent in the perpendicular and parallel 
directions known as "critical balance" may be derived. 
To understand this feature, consider the timescale T\, 
the energy-transfer time at the corresponding scale A, 
characterizing the nonlinear dynamics at that scale. T\ 
does not necessarily coincide with the eddy turnover time 
t\ = X/Szx because of the Alfven effect. Spatial struc- 
tures along the axial direction result from wave propaga- 
tion (at the Alfven speed ca) of the orthogonal fluctua- 
tions. In other words, the cascading of turbulence in two 
different planes separated by a distance In leads to for- 
mation of sc ales in the parallel direction whose smallest 
size can be (iGoldreich fc Sridharf Il995t ICho efail 120021 : 
lOughton et al.lll994h 



*||(A) ~c A T x , 



(44) 



the critical balance condition. T\ will be smaller at 
smaller scales, so that smaller perpendicular scales create 
smaller axial scales. 

Figure [9] shows a snapshot at time t ~ 145 T4 of the 
2D spectrum E(n±,n z ) for run I in bilogarithmic scale, 
where n z and n± are respectively the axial and orthogo- 
nal wavenumbers. Consider vertical cuts at n± = const: 
it is clearly visible that from n± = 1 up to n± ~ 20 the 
wavenumbers with n z > 1 are scarcely populated com- 
pared to the respective wavenumbers with n z < 1 (the 
parallel spectrum has also the n z = component, in Fig- 
ure[9]the vertical coordinate is n z + l). However note also 
how the loci of maximum parallel wave-number do not 
precisely follow the critical balance line, rather they are 
offset at larger n±: in our case, the hypothetical length of 
the axial structures (from critical balance) can be longer 
than the characteristic length of the system, in our case 
the length of the coronal loop L. But in the range of 
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Fig. 10. — : Total energy at the injection scale (modes 
3 and 4), time-averaged for the four simulations F, G, H 
and I with different Alfven velocities. The dashed line 



shows the curve E*, 



cx c 



A- 



while the continuous line 



shows Ei n as a function of c A as obtained from equa- 
tion (|66[) for a = corresponding to a Kolmogorov spec- 
trum. The actual growth of £ m , both simulated and 
derived from (|66|) . show that the growth is less than 
quadratical but higher than in the simple Kolmogorov 
case. 



perpendicular wavenumbers for which 

*I|(A)>L, (45) 

boundary conditions, i.e. line-tying, intervenes and the 
cascade along the axial direction is strongly inhibited. In 
our simulations this occurs roughly at n± ~ 20. Beyond 
n± ~ 20 the spectrum is roughly constant along n± — 
const up to a critical value where it drops. 

Interestingly enough, the slope of the ID spectrum for 
run I (Figure [8]) diminishes its value around n± ~ 20. 
The reason is that the condition ^ii(A) > L with £i\ (A) 
defined by critical balance, turns out to play a major 
role in the "strength" or "weakness" of the cascade: for 
n±_ < 20 the system is in a weak turbulent regime, while 
for n± > 20 a transition to strong turbulence is observed. 

In our runs, larger values of c^t, i.e. of the parameter 
/ (fT8|) . lead to larger magnetic energy and total energies, 
while the kinetic energy remains smaller than magnetic 
energy and increases much less (increasing its value by a 
factor of 6 from c_\ = 50 to c_4 = 1000). In particular 
Figure flOl shows total energy at the injection scales (see 
§ 12. 2p . i.e. the sum of the modes n = 3 and n = 4 (see 
eq. |43| of total energy, 

E m = E 3 + E 4 (46) 

as a function of the non-dimensional Alfven velocity c^. 
Their growth is less than quadratical in c^, which implies 
that the rms of the velocity and magnetic fields at the 
injection scale (or equivalently the Elsasser fields Szi n ) 
grow less than linearly. Hence as c A is increased, the ra- 
tio x, a measure of the relative strength of the nonlinear 
interactions at the injection scale, decreases: at differ- 
ent values of different regimes of weak turbulence are 
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therefore realized at the larger scales of the inertial range, 
as the different spectra in Figure [8] confirm. 

The presence of a "double" inertial range, with a weak- 
type power-law index at larger scales, and a flatter strong- 
type power-law index at smaller scales would not affect 
the overall cascade rate, and therefore the scalings of 
loop heating with loop parameters. These are set at the 
larger scales, and are therefore dependent on the cascade 
rate determined by the weak-type scaling law, for which a 
physically motivated phenomenological derivation is pre- 
sented below. We stress that the possible existence of a 
"double" inertial range, surmised here with scaling laws 
and somewhat tenuous numerical evidence, does not ap- 
pear to have been predicted before and requires substan- 
tiating evidence from higher numerical resolution simu- 
lations which are planned for the near future. 

5.2. Phenomenology of the Inertial Range and Coronal 
Heating Scalings 

We now introduce a phenomenological model to deter- 
mine the energy transfer time-scale T\ and as a conse- 
quence the properties of the cascade. This time-scale, 
and therefore the different spectra which result, can only 
depend on the single non-dimensional quantity defining 
our system, namely / = £ c vaI Lu p h (|18p . The simula- 
tions show that as this parameter is increased the spectra 
steepen leading to a weakened cascade. We revert here 
to dimensional quantities for the scalings, so that we can 
quantify the resulting coronal heating rates. 

The Alfven effect is based on the idea that two coun- 
terpropagating Alfven waves interact only for the time 
tii = £\\/va, leading to a transfer energy time longer that 
the "generalized" eddy turnover time 



A 

5z\ ' 



(47) 



X 



(48) 



The ratio between these two timescales 

T A _ t\\ 5z x 

ives a me asure of th e ir rel ative strength. Ilroshnikovl 
1964) and iKraichnanl (|1965| ) proposed that the energy 
transfer time T\, because of the Alfven effect, is longer 
than the eddy turnover time, and is given by 



rp t a 



(49) 



where however they considered an isotropic situation, so 
that the Alfven time was given by the propagation time 
over the scale of the Alfvenic packet. For weak turbu- 
lence however i\\ > L, so that the Alfven time must be 
based on the scale L: ta = L/v A . 

In addition, we must allow line-tying which acts to slow 
the destruction of eddies on a given scale A more effec- 
tively than the standard random encounter effect t\ / ta 
(jDobrowolnv et ail Il980h . We can therefore assume a 
sub-diffusive behaviour for z + z non-linear en- 
counters leading to 



J A 



(50) 



with values a > 1 and depending in some way on the pa- 
rameter / [recall that a = 0, 1 correspond respectively to 



anisotropic Kolmogorov and Kraichnan cases (the latter 
leading to a fc~ 2 inertial range spectrum)]. 

Our simulations then close this ansatz by determining 
how alpha depends on /: integrating over the whole vol- 
ume (£ x £ x L), the energy cascade rate may now be 
written as 

x 2 



Lp- 



Using (|50|) the energy transfer rate is given by 



L ■ p — ^ 



l> \ — 

VA 



v Q+3 



A Q + ! 



(51) 



(52) 



Identifying, as usual, the eddy energy with the band- 
integrated Fourier spectrum Sz 2 ~ k^Ej t± , where k± ~ 
£/X, from eq. (|52[) we obtain the spectrum 

E k± <xkJ^, (53) 

where for a = 0,1 the —5/3,-3/2 slope for the 
anisotropic Kolmogorov, Kraichnan spectra are is recov- 
ered, but steeper spectral slopes up to an asymptotic 
value of —3 arc obtained with higher values of a. 

Correspondingly, from eqs. (|5Tj) - (f5"2"|) , the following 
scaling relations for Szx and Ta follow: 



Sz x 



\PLp 
£ 2 Lp 



L 



i-3 



L 



1+3 . 2+1 

A=+ 3 



A 



(54) 
(55) 



Recently iBoldvrevI ([2005) has proposed a similar 
model, which aims to overcome some discrepancies be- 
tween previous models and numerical simulations, and 
that self-consistently accounts for the formation of cur- 
rent sheets, for the cascade of strong turbulence. His 
energy transfer time is given by 

A / va \ a 



Tx = 



Sz x \SzxJ 



(56) 



but he suggests the interval < a < 1 as appropriate to 
strong turbulence. 

As pointed out above of § [3J the solutions of equa- 
tions (fT2 |) -(fT4 | depend only on the non-dimensional pa- 
rameter / = £ c va/ Lu p h (eq. (fT8]) ) and so a (|50|) is only 
a function of / 

LVA ^ (57) 



a 



Lu p h 



We estimate the value of a from the slope of the to- 
tal en ergy spectra (l53|) . as described in iRappazzo et all 
(|2007h . As shown in Figure [8] to different values of 
CA = VA/uph, (i.e. /) ranging from 50 up to 1000 corre- 
spond spectral slopes from ~ —1.8 up to ~ —2.7. Thes in 
turn correspond (through cq. (|53p ) to values of a ranging 
from ~ 0.33 up to ~ 10.33. 

How do the above results affect coronal heating scal- 
ings? The energy that is injected at the large scales 
by photospheric motions, and whose energy rate (em) is 
given quantitatively by the Poynting flux (|2ip . is trans- 
ported (without being dissipated) along the inertial range 
at the rate e (l52)l . to be finally dissipated at the rate €d- 
In a stationary state all these fluxes must be equal 

tin = e = e d (58) 
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The injection energy rate (|21[) is given by S, the Poynt- 
ing flux integrated over the photospheric surfaces: 



e in — S — 
= PVA 



[ da {u\ ■ b±) - f da (u° ± ■ b ± ) 

Jz=L Jz=0 



(59) 



2D spatial periodicity in the orthogonal planes allows 
us to expand the velocity and magnetic fields in Fourier 
series, e.g. 



U-L {x,y) = ^u r , s 



where 



2-7T 

K,s = -j- (r, s, 0) r,seZ 



(60) 



(61) 



The surface integrated scalar product of uj^ and at 
the boundary is then given by 

Jdau^ ■ b± = ^2 u r\s ■ J J dxdy bj_e lkr - s ' x = 



, s e Z (62) 



This integral is clearly dominated by large scales con- 
sistent with observations of photospheric motions. In 
our case (eq. (|10j0 boundary velocities only have compo- 
nents for wave numbers (r, s) € Z 2 with absolute values 

between 3 and 4, 3 < (r 2 + s 2 ) 1/2 < 4. Then in §2§ 
only the corresponding components of b^ are selected. 

At the injection scale, which is the scale of convective 
motions £ c ~ 1,000 fern, a weak turbulence regime de- 
velops, so that the cascade along the axial direction z is 
limited and in particular the magnetic field b^ can be 
considered approximately uniform along z at the large 
orthogonal scales. Then from eq. (fB"9"|) we obtain 



= S 



PVA 



da (u^ 



b, 



(63) 



Introducing u p h = — Uj_, using relation ([62]) . and 
integrating over the surface, we can now write 



S ~ l 2 pvAU ph 5ze c 



(64) 



where we have approximated the value of Sbi c , the rms 
of the magnetic field at the injection scale £ c , with the 
rms of the Elsasser variable because the system is mag- 

1/2 

netically dominated, i.e. 5zg c — [5u\ + 5b 2 ) ~ Sbg c . 

We now have an expression for t in , where the only 
unknown variable is 6zg c , as £ c , p, va and u p h are the 
parameters characterizing our model of a coronal loop. 

The transfer energy rate e does not depend on A. Con- 
sidering then A = i c in equation |52|) . we have 



pl 2 L a+1 



v a+3 



ea+l , 



(65) 



Equations ([64]) and ([65]) show another aspect of self- 
organization. Both ti n and e, respectively the rate of the 
energy flowing in the system at the large scales, and the 
rate of the energy flowing from the large scales toward 
the small scales depend on Sz£ c , the rms of the fields at 
the large scale. This shows that the energetic balance 



of the system is determined by the balance of the en- 
ergy fluxes e and ei„ at the large scales. The small scales 
will then dissipate the energy that is transported along 
the inertial range (see eq. (|58|) ). This implies that be- 
yond a numerical threshold total dissipation (dissipation 
integrated over the whole volume) is independent of the 
Reynolds number. In fact beyond a value of the Reynolds 
number for which the diffusive time at the large scale is 
negligible, i.e. when the resolution is high enough to re- 
solve an inertial range, the large-scale balance between 
e and ej„ is no longer influenced by diffusive processes. 
Of course this threshold is quite low respect to the high 
values of the Reynolds numbers for the solar corona, but 
it is still computationally very demanding. 

An analytical expression for the coronal heating scal- 
ings may be obtained from (|64[) and (|65p yielding the 
value of bz* for which the balance £,„ = e is realized: 



\ Lup h J 



»+2 



(66) 



Uph 

Substituting this value in (|65p or equivalently in (1641) we 
obtain the energy flux 



S* 



( icVA \ 



ipVAU ? h [tu-J 



(67) 



As stated in (|58|) in a stationary cascade all energy fluxes 
are equal on the average. S* is then the energy that for 
unitary time flows through the boundaries in the coronal 
loop at the convection cell scale, and that from these 
scales flows towards the small scales. This is also the 
dissipation rate, and hence the coronal heating scaling, 
i.e. the energy which is dissipated in the whole volume 
for unitary time. As shown in equation (|57p the power a 
depends on the parameters of the coronal loop, and its 
value is determined numerically with the aforementioned 
technique. 

The observational constraint with which to compare 
our results is the energy flux sustaining an active region. 
The energy flux at the boundary is the axial component 
of the Poynting vector S z (see § 13. ip . This is obtained 
dividing S* (|6"7) , the Poynting flux integrated over the 
surface, by the surface i 2 : 



z = ~p ~ P v AU ph 



£ C VA 



Lu 



ph 



1 + 2 



(68) 



where a is not a constant, but a function of the loop 
parameters (|5T|) . The exponent in (f6"B")) goes from 0.5 for 
a = up to the asymptotic value 1 for larger a. We 
determine a numerically, measuring the slope of the in- 
ertial range (Figure [5]) , and inverting the spectral power 
index (|53p . We have used simulations F, G, H and I to 
compute the values of a, because they implement hy- 
perdiffusion, resolve the inertial range, and then are be- 
yond the numerical threshold below which total dissipa- 
tion does not depend on the Reynolds number. These 
simulations implement va — 50, 200, 400 and 1,000, and 
the corresponding a are ~ 0.33, 1, 3, 10.33. The corre- 
sponding values for the power (a + l)/(a+ 2) (f6"8")) are 
~ 0.58, 0.67, 0.8 and 0.91, close to the asymptotic value 
1. S z is shown in Figure fTTI (diamond points) as a func- 
tion of the axial Alfven velocity va- To compute the 
value of S z for va = 2, 000 fcms -1 we have estimated 
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Fig. 11. — : Analytical ()68|l and numerically computed 
dissipated flux as a function of the axial Alfven velocity 
Va- The continuous line shows the Poynting flux (|68[) 
as a function of va in the case a = 0, corresponding 
to a Kolmogorov-like cascade. To higher values of va 
correspond a higher dissipation rate because a weak tur- 
bulence regime develops. 



a ~ .95, although for values close to 1 S z does not have 
a critical dependence on the value of the exponent. 

In Figure[Tl]we compare the analytical function S z (|68|) 
with the respective value determined from our numer- 
ical simulations (star points), i.e. with the total dis- 
sipation rate by the surface and converted to dimen- 
sional units ((J + Vl)/l 2 , see (|37)l ). For the numerical 
simulation values, the error-bar is defined as 1 stan- 
dard deviation of the temporal signal. The analytical 
and computational values are in good agreement for all 
the 4 simulations considered, and for the more realis- 
tical value va = 2, 000 ferns -1 the dissipated flux is 
~ 1.6 x 10 6 erg cm s -1 . This value is in the lower range 
of the observed constraint 10 7 erg cmT 2 s -1 . 

The continuous line in Figure [IT] corresponds to the 
function S z for a = (which is approximately realized 
for va < 50 km s -1 ), in correspondence of which a Kol- 

3/2 

mogorov spectrum would be present, and S z oc vl . The 
computed and analytical values of S z for higher va are 
always beyond this curve, because a increases its values, 
and a more efficient dissipation takes place. This is due 
to the fact that to higher values of a correspond higher 
values of the energy transfer time, and consequently a 
longer linear stage, higher values of the fields at the large 
scales (|66|) , and hence a higher value of the energy rates 
(see (|64"1) . ([65]) and (|67|) ). So that it is realized the only 
apparently paradox that to a weaker turbulent regime, to 
which corresponds less efficiency in the nonlinear terms, 
corresponds a higher total dissipation. 

In the last paragraph of § l3.1l we have shown that when 
the condition (|2"3")l is satisfied the emerging flux can be 
neglected. But in eq. (|23|) we have to specify the value of 
the magnetic field b l ^ rb self-consistently generated by the 
non- linear dynamics. This value is given by (|66p as the 
magnetic field dominates (Sz^ ~ By substitution 

we can now estimate that the emerging flux is negligi- 
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Fig. 12. — : Transition to turbulence: Total ohmic 
and viscous dissipation as a function of time for simula- 
tions A, B, C, G (displayed on the same scales). All the 
simulations implement ca = 200, but different Reynolds 
numbers, from Re = 200 up to 800. Run G implements 
hyperdiffusion. For Reynolds numbers lower than 100 
the signal is completely flat and displays no dynamics, at 
higher Reynolds smaller temporal structures are present. 

ble when the emerging component of the magnetic field 
satisfies 



(69) 



bi < Bailie 




In the asymptotic state a » 1 the condition reduces 
to b*f / Bo < y/£ c /L. For a coronal loop with L ~ 
40, 000 km, as i c ~ 1,000 fern this implies that emerg- 
ing flux does not play a role if b e / / Bq < 1/6. 

5.3. Transition to Turbulence and Dissipation vs. 
Reynolds Number 

Turbulence is a ch a racter istic of high Reynolds number 
systems (e.g.. iFrischl (|1995l) ). For a sufficiently high vis- 
cosity nonlinear dynamics is strongly suppressed, and our 
system relaxes to a diffusive equilibrium ( £I3.3|) . and no 
significant small scale is formed. Increasing the Reynolds 
number, the diffusive time at the injection scale (JHJ) 
Td r~j Re l 2 c increases. At a certain point it will be big 
enough not to influence the dynamics as the large scales, 
an inertial range will then be resolved and total dissipa- 
tion will not depend any longer from the Reynolds num- 
bers. In fact for higher values of Re the intertial range 
will extend to higher wave-numbers, but the energy flux 
will remain the same. 

At higher Reynolds numbers smaller scales are re- 
solved, and each scale will contribute with its charac- 
teristic time T\ to the temporal structure of the rms 
of the system. Figure [12] shows total dissipation as a 
function of time for simulations A, B, C and G, on the 
same time interval, and on the same scale. At increas- 
ingly higher values of the Reynolds numbers smaller and 
smaller temporal structures are added to the signal. Ide- 
ally the temporal structure of total dissipation at higher 
Reynolds numbers is well described by shell-model simu- 
lations. For smaller values of Re the signal is completely 
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Fig. 13. — : Total ohmic and viscous dissipation as a 
function of time for simulations A, B, C, D, E and G, all 
of them implement ca = 200 but different Reynolds num- 
bers. The threshold beyond which dissipation is indepen- 
dent of the Reynolds number can be identified around 
Re = 800, corresponding to a numerical resolution of 
512x512 points in the orthogonal planes. 

flat (see Figure [T3]) . This behaviour identifies a transition 
to turbulence. 

Figure [13] shows total dissipation as a function of time 
for the same 4 simulations shown in Figure [T2l plus other 
2 simulations at lower Reynolds number, respectively 
Be = 100 and Re — 10 for the complete time interval. 
For the lowest value of Re no dynamics is present, so 
that the threshold value for the transition to turbulence 
can be set to Re ~ 100. To higher values of Re dissi- 
pation grows. An inertial range is barely solved with a 
resolution of 512x512 grid points in the x-y planes, so 
that simulation with Re = 800 can be considered at the 
threshold. On the other hand simulation G implements 
hyperdiffusion, so that an inertial range is solved, and 
the dynamics is not affected by diffusion. The presence 
of a sufficiently extended inertial range implies in fact 
that we are beyond the numerical threshold where dissi- 
pation does not depend on the Reynolds number (§ 15. 2p . 
The threshold value can be identified to a sufficient ex- 
tent at Re = 800, i.e. for a numerical grid of 512x512 
points. The number of points to use along the axial di- 
rection should be enough to allow the formation of all the 
small scales due to the "critical balance" (Figure O, but 
a larger number of points would only result in a waste of 
computational time. 

5.4. Inverse Cascade and Line-tying 

Two dimensional simulations (|Einaudi et alj 119961 : 
iGeorgoulis et al.lll998[ ) have shown an inverse cascade for 
the magnetic energy, corresponding in physical space to 
the coalescence of magnetic islands. In the 3D case the 
DC magnetic field along the axial direction is present, 
giving rise to a field line tension that tends to inhibit 
an inverse cascade, as motions linked to the coalescence 
would bend the field lines of the total magnetic field, 
which are mostly elongated along the axial direction (Fig- 
ure 18). On the other hand field line tension depends on 
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Fig. 14. — : Run F: In this simulation with ca = 50 
an inverse cascade at the wavenumbers n = 1 and 2 is 
realized. Energy is injected at wavenumbers n = 3 and 
4. 



the strength of the axial field, becoming stronger for a 
stronger field. 

In Figures Q]|] and [15] the first 4 modes of magnetic en- 
ergy for simulations F and I, respectively with ca = 50 
and 1000, are plotted as a function of time. Energy is 
injected at wave-numbers n ~ 3 and 4. Modes associated 
to wave-numbers 1 and 2 grow to higher values than at 
the injection scale in run F, while in run I they are always 
limited to lower values. In runs G and H, with respec- 
tively ca — 200 and 400 an intermediate behaviour is 
found, but none of the modes n = 1 or 2 never becomes 
bigger than the injection energy modes. 

5.5. Timescales 

In the previous sections we have always affirmed that 
the Alfven crossing time ta = L/va is the fastest 
timescale in the system, and that in particular it is 
smaller than the nonlinear timescale T n i, which we can 
identify with the energy transfer time (|55[) at the injec- 
tion scale T n i — Ti c . 

In Figure [3] it is already clear that the nonlinear 
timescale is longer that ta, in fact it shows that the 
timescale over which energy has substantial variations is 
bigger than the Alfven crossing time. 

The same behaviour is identified in Figures fTHTrSl 
which show the time evolution of the magnetic energy 
modes for runs F and I. These are more relevant quanti- 
ties, because to realize a weak MHD turbulence regime 
it is required that the energy transfer time T\ is bigger 
than the crossing time ta at the injection scale A = £ c 
and for a limited range of smaller scales down to some 
lower bound A*: A* < A < £ c . The magnetic energy 
modes at the injection scale (n = 3 and 4) change their 
values on scales bigger than ta, and for a larger value 
of the Alfven velocity the nonlinear timescale is longer 
respect to the crossing time (Figures [T4l[T5|) . We can 
roughly estimate r„; ~ 5ta for run F with ca — 50 and 
r n i ~ 20ta for run I with ca = 1, 000. 

Using our scaling relations we can derive an analytical 
estimate for the energy transfer time T\. Substituting 
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Fig. 15. — : Run I: Simulation performed with C4 = 
1000. The increased magnetic field line tension inhibits 
an inverse cascade for the orthogonal magnetic field. 



the energy rate (|67|) in equation (|55j) we obtain: 

Q Q + l 



ir 



(70) 



where r c = £ c /u p h- In particular the ratio over the Alfven 
crossing time is: 



Tx 

TA 



1+1 

i+2 



(71) 



and as t c > T4 then self-consistently T\ > T4. For 
our loop £ c ~ l,000fcm and itp/j ~ lfcms^ 1 , so that 
r c ~ 1,000 s. For runs F and I shown in Figures [T4l and 
[TBI the loop length is always L = 40, 000 km, while the 
Alfven velocity is respectively va — 50 and 1, 000 km a , 
and the corresponding crossing times T4 = 800 and 40 s. 
Using the values of a computed in £15.21 (respectively 
a = 0.33 and 10.33) we can then roughly estimate from 
(|7ip . the nonlinear timescale r„; = T\ = i c and its ratio 
with the Alfven crossing time: 



1~nl 
TA 



TA 



1 + 1 



(72) 



For runs F and I we find t„;/t4 = 1.2 and 22.3 in agree- 
ment with the simulations. 

Equation (171]) can also be used to estimate the exten- 
sion of the weak turbulence inertial range. The region for 
which the weak turbulence condition T\ > ta is satisfied 



is: 



\> A* =£ r \ — 
T r 



(73) 



Figure [TB] shows the temporal spectrum of magnetic 
energy for run G with ca = 200, i.e. we perform the 
Fourier transform of the magnetic energy as a function 
of time, and then plot its squared modulus. We use 
run G because it is the one for which we have saved 
more frequently the rms quantity and then the plot cov- 
ers a wider range at high frequencies. The power spec- 
trum is roughly constant up to vJva ~ 0.2, which cor- 
responds to t/rA ~ 5 in agreement with our scaling (|T2[) 
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Fig. 16. — : Temporal spectrum of magnetic energy for 
run G. va — 1/ta is the frequency corresponding to the 
Alfven crossing time. The intermediate part of the spec- 
trum exhibits a v~ 2 power law. 



which for this case gives t„;/t4 ~ 3.3. Beyond this 
critical point the power spectrum exhibits a power law 
which fits v~ 2 , in agree ment with shell-model simulations 
ijBuchlin fc Vellil[200l . 

6. DISCUSSION AND CONCLUSIONS 

We would like first to clarify a few concepts that might 
otherwise result in misunderstandings of the work that 
we have presented. The concept of turbulence is used 
to describe different processes in different research fields, 
so that its use, without specifications, can result vague 
and misleading. It is in fact very often used to describe 
chaotic behaviors at the small scales, often linked to the 
intermittent dissipation of energy. Although this aspect 
is present in our simulations, when we say that the Parker 
problem is an MHD turbulence problem, we refer mainly 
to the property of turbulence to transfer energy from 
large to small scales. Namely to its ability to trans- 
port the energy from the scale of photospheric motions 
(~ 1000 km), where it is injected, down to the small dis- 
sipative scales (meters?), without dissipating it at the 
intermediate scales. This property is clearly identified 
by the presence of an inertial range with a power law 
spectrum, which extends from the injection scale to the 
dissipative scale. 

Furthermore turbulence, magnetic reconnection and 
ohmic heating associated to currents are sometime pre- 
sented as alternative and/or mutually exclusive coronal 
heating models. This contraposition is artificial. Cur- 
rent sheets are in fact the dissipative structures of 
MHD turbulence, and magnetic reconnection at the loci 
of current sheets is observed in virtually every MHD tur- 
bulenc e simulation in both 2D and 3D (see e.g.. lBiskamd 
(2003) and references therein ). Nanoflares are then nat- 
urally associated with the time and space intermittency 
of the sma ll scale deposition of en ergy (as shown in the 
2D case bv iGeorgoulis et ail (|1998h ). which is due to the 
cascade which leads to the formation and dissipation of 
current sheets, and to which we refer collectively with 
the term MHD turbulence. 
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In summary, the main results presented in this paper 
are the following: 

• The time-dependent Parker problem may be seen 
as an MHD turbulence problem, where the large 
scale forcing function is realized by the photo- 
spheric motions. 

• This system is genuinely turbulent, in the sense 
that small scale formation is not driven passively 
by the random walk of the footpoints, rather it 
is a property of the maxwell stresses developing 
in the coronal volume. Current sheets therefore 
do not generally result directly from a "geomet- 
rical" misalignment of neighboring magnetic field 
lines stirred by their footpoint (random) motions, 
they are the result of a nonlinear cascade in a self- 
organized system. 

• Nanoflares are naturally associated with the inter- 
mittent dissipation of the energy that, injected at 
the large scales by photospheric motions, is trans- 
ported to the dissipative scales through a cascade, 
and is finally dissipated through nonlinear mag- 
netic reconnection. 

• Beyond a threshold, that is low compared to the 
coronal Reynolds numbers, but still computation- 
ally very demanding, total dissipation is indepen- 
dent of the Reynolds numbers. This threshold cor- 
responds to a numerical resolution of ~ 512 x 512 
grid points in the planes orthogonal to the domi- 
nant DC magnetic field. 

• As the loop parameters vary, different regimes of 
turbulence develop: strong turbulence is found for 
weak axial magnetic fields and long loops, leading 
to Kolmogorov-like spectra in the perpendicular di- 
rection, while weaker and weaker regimes (steeper 
spectral slopes of total energy) are found for strong 
axial magnetic fields and short loops. There is no 
single universal scaling law (see (|68jl ). as a con- 
sequence the scaling of the heating rate with ax- 
ial magnetic field intensity, which depends on the 
spectral index of total energy for given loop param- 
eters, must vary from B^ 2 for weak fields to B 2 for 
strong fields at a given aspect ratio. 

• For a loop 40, 000 km long , with an Alfven veloc- 
ity Va = 2,000 km s" 1 and a numerical density of 
10 10 cm" 3 , whose footpoints are subject to photo- 
spheric motions of u p h ~ lkms" 1 on a scale of 
l c ~ 1, 000 km, the energy flux entering the system 
and being dissipated is S z ~ 1.6 x 10 6 erg cm" 2 s" 1 . 
On the other hand, for a coronal loop typical of a 
quiet-Sun region, that has the same parameter of 
the previous case but with a length of 100, 000 km 
and va — 500 kms , the resulting Poynting flux 
is S z ~ 7 x 10 4 erg cm" 2 s" 1 . 

The most advanced EUV and X-RAY imagers (e.g. 
those onboard SOHO, TRACE, STEREO and HINODE) 
have space resolutions (~ 800 km) of the order of the 
granulation cells. Hence they do not resolve the small- 
scales where current sheets, magnetic reconnection and 
all the dynamical features of the system take place. Their 



resolution is roughly 1/5 the length of the perpendicular 
cross-section of our numerical box (~ 4000 km) . Hence, 
even if the system is highly dynamical on small-scales 
(see Figure 17 and the associated movie), integrating 
over these scales has the effect to "averaging" the small 
scale dynamics. In particular small scale reconnection 
cannot be detected, magnetic fieldlines will appear only 
slightly bended (Figure 18), and their dynamics will ap- 
pear slower (a modulation of the nonlinear timescale with 
the thermodinamical timescales). 

The topological and dynamical effects associated with 
magnetic reconnection should be taken into account 
when modeling the thermody namical and ob servational 
properties of coronal loops (jSchriiverl l2007f) . recalling 
that most of the dynamics take place at sub-resolution 
scales while we observe the integrated emission. 

Two density current fields that have the same "steady" 
integrated ohmic dissipation, balanced by a correspond- 
ing Poynting Flux (see § 13.31 equation (|37|) and Figured]), 
but with different spatial distributions will have different 
emissions. Consider the first with only large scale com- 
ponents, as the one that would result from a diffusive 
process (§ 13. 3|) . while in the second the current has only 
small scale components, as in the simulations that we 
have presented. In the second case the filling factor is 
small (Figures 17 and 18) so that the density of current 
has a far larger value, and this would correspond to two 
very different thermodynamical and observational out- 
comes. But the highly dynamical effects associated with 
the second case will be averaged and result less dynamical 
when integrated. Still the integrated observables should 
be very distinct between the two cases. 

Finally, while our simulations give an accurate descrip- 
tion of the time-dependent Parker problem, with the lim- 
itations on the photospheric forcing field described in the 
introduction, the use of the reduced MHD equation is 
justified only for slender loops threaded by a strong ax- 
ial magnetic field. For short loops, or loops that have 
orthogonal component of the magnetic field compara- 
ble to the axial component, the full set of MHD equa- 
tion should be implemented. For the slender loops that 
we have simulated we observe a modest accumulation 
of energy, which subsequently is released via nanoflares. 
On the other hand shorter loops, or loops in a more 
complicated geometry, or subject to loop-loop interac- 
tions, and more generally loops affected by the neigh- 
boring coronal environment, might exhibi t the ability to 
accumulate more energy (e.g. Low (2006)) and then re- 
lease itJ^J^jge£flBres_ 2 _gossibly via a "secondary instabil- 
ity" dD ahlbur g et al.ll2005h or fast magnetic reconnection 
()Cassak et al.ll2006t ). 

Magnetohydrodynamics (MHD) has proved to be a 
useful to ol to investigate the properties of the turbulent 
cascade (|Biskampl l200l . MHD is very well known to 
give an approximate description of the plasma dynamics 
at large scales and low frequencies. In MHD turbulence 
it is generally supposed that at the small scales a "dis- 
sipative mechanism" is present. Most of the properties 
of the turbulent cascade do not depend on the details of 
the dissipative mechanism, whether it is described by the 
diffusive operator present in equations ([T])-([2|), or more 
properly by a kinetic mechanism. 

In particular in our case, the timescales associated at 
the scale A ( ([70]) for weak turbulence and (|56|) for the 
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strong case) decrease for smaller scales. In this way the 
small-scale dynamics is characterized by high-frequency 
phenomena, and then it is not well described by MHD, 
but rather a kinetic model would be more appropriate. It 
is then possible that (self-consistently) at the small scales 
particle acceleration plays an important role in the dissi- 
pation of energy, a physical process that should be inves- 
tigated through kinetic models. Nevertheless the coronal 
heating rates , like the cascade properties over an ex- 
tended range of scales, are independent of the details of 
the dissipation mechanism. They are determined by the 
balance, at the large scales (see § I5.2p . between the rate 
of the energy flowing into the loop from the boundaries 
due to the work done by photosphcric motions on the 
magnetic field line footpoints at the scale of the convec- 
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tive cells, and the rate at which the energy flows along 
the inertial range from the large scales towards the small 
scales. 
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Fig. 17. — : Run A: (a) Streamlines of the boundary velocity fields u 1 ^ — constant in time, (b)-(e) Axial component 
of the current j (in color) and field- lines of the orthogonal magnetic field in the mid-plane (z = 5), at selected times 
covering the linear and nonlinear regimes up to t = 18.47 T4. (f) Axial component of the vorticity lu (in color) and 
field-lines of the orthogonal magnetic field in the mid-plane at time t — 18.47 r^. 

During the linear stage the orthogonal magnetic field is a mapping of the boundary forcing (cfr. a and b). After the 
collaps of the large-scale currents (b, c, d), which in Fourier space correspond to a cascade of energy (see Figure [5]), 
the topology of the magnetic field departs from the boundary velocity mapping and evolves dynamically in time 
(see movie), (e)-(f) Current sheets are embedded in quadrupolar vorticity structure, a clear indication of nonlinear 
magnetic reconnection. 
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Fig. 18. — : Run A: Side (a,c) and top (b,d) views of cu rren t sheets (a,b) and field lines of the total magnetic field 
(c,d) at time t = 18.47 T4 (same time as in Figures fl7clll711) . For an improved visualization the box size has been 
rescaled, but the axial length of the computational box is 10 times longer that the perpendicular cross-section length. 
The rescaling of the box artificially enhances the structures inclination. To restore the original aspect ratio the box 
should be stretched 10 time along z. 

(a)-(b) Two isosurfaces of the squared current j 2 . The isosurface at the value j 2 = 2.8 • 10 5 is represented in partially 
transparent yellow, while red displays the isosurface with j 2 — 8 • 10 5 , well below the maximum value of the current 
at this time j^ax = 8.4 • 10 6 . As typical of current sheets, isosurfaces corresponding to higher values of j 2 are nested 
inside those corresponding to lower values. For this reason the red isosurface appears pink. Although from the side 
view the sheets appear space-filling, the top view shows that the filling factor is small. 

(c)-(d) Field-lines of the total magnetic field (orthogonal plus axial), and in the mid-plane (z = 5) field- lines of the 
orthogonal component of the magnetic field. 



